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Abstract 



O 

PsJ ' We investigate the fuU functional form of the first passage time density (FPTD) of a tracer particle 

in a single-file diffusion (SFD) system whose population is: (i) homogeneous, i.e. all particles having 



the same diffusion constant and (ii) heterogeneous, with diffusion constants drawn from a heavy- 
tailed power-law distribution. In parallel, the full FPTD for fractional Brownian motion [fBm - 



■ defined by the Hurst parameter, H G (0, 1)] is studied, of interest here as fBm and SFD systems 

Q . belong to the same universality class. 

• 1— I ■ 

\ Extensive stochastic (non-Markovian) SFD and fBm simulations are performed and compared 

' to two analytical Markovian techniques: the Method of Images approximation (MIA) and the 

*^ \ 

Willemski-Fixman approximation (WFA). We find that the MIA cannot approximate well any 

^ ■ 

, temporal scale of the SFD FPTD. Our exact inversion of the Willemski-Fixman integral equation 

captures the long-time power-law exponent, when H > 1/3, as predicted by Molchan [1999] for fBm. 
When H < 1/3, which includes homogeneous SFD [H = 1/4), and heterogeneous SFD [H < 1/4), 
On . the WFA fails to agree with any temporal scale of the simulations and Molchan's long-time result. 

V/^ ' SFD systems are compared to their fBm counter parts; and in the homogeneous system both 

o ■ 

CN ' scaled FPTDs agree on all temporal scales including also, the result by Molchan, thus affirming 

that SFD and fBm dynamics belong to the same universality class. In the heterogeneous case 



^ ! SFD and fBm results for heterogeneity- averaged FPTDs agree in the asymptotic time limit. The 

non-averaged heterogeneous SFD systems display a lack of self-averaging. An exponential with a 
power-law argument, multiplied by a power-law pre- factor is shown to describe well the FPTD for 
all times for homogeneous SFD and sub-diffusive fBm systems. 



* Corresponding author: |lloyd.sanders@thep.lu.se| 
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I. INTRODUCTION 



Within the physical phenomena that exhibit stochasticity, the idea of the first passage 
time density (FPTD) is of great importance Recently the literature has seen analysis 
of FPTDs in a wide variety of research areas, including stochastic systems: Diffusion with 
stochastic resetting diffusion in complex scale-invariant media cosmology (the cal- 
culation of the mass distribution of dark matter halos) [4]; evolutionary ecology (foraging 
tactics of various animal species) jsl; and covered extensively in primarily one-dimensional 
(ID) single particle scenarios by Redner [1]. Most previous studies focus on FPTD for 
Markovian (memory-less) dynamics. 

The role of crowding in physical systems is of key relevance in providing non-Markovian 



An important special case 



effects in many scenarios, especially within the cell (e.g. 
in this area of stochastic processes is single-file diffusion (SFD), defined as many particles 
diffusing in a ID system, where the order of the particles within the system remains constant, 
i.e. the particles are subject to hard-core repulsion. Experimentally SFD is shown to exist in 
many types of systems, such as the diffusion of colloids in channels 0, lo|, and of relevance 
in the motion of fiuorescently tagged proteins on DNA 

Initially, because of the problem's mathematical tractability; over the past five decades, 
beginning with the hallmark theoretical investigations of Harris [12;], followed by Lebowitz 
and Percus fisl, and, for instance. Fisher jl4|, a large amount of (theoretical) work has been 
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20|, generalizing the original 



conducted into SFD systems, notably in the latter decade 
results to particles interacting with general potentials and in external force fields. Of interest 
for the present study is the recent increasin g a mount of research into complex population 



types or heterogeneous SFD systems 
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2ii-y|. 



One prominent example where SFD and the FPTD are amalgamated is within the context 
of protein-DNA interaction. It was shown that the time taken for a protein to find a specific 
start location is up to two orders of magnitude greater than predicted by Smoluchowski 
three-dimensional (3D) diffusion rate alone j25|. This led to the concept of "facilitated" 



target location - the combination of 3D diffusion through the cytop 
along the DNA; this coupling is shown to exist experimentally 



theoretically 



26 



asm and ID diffusion 



28|, and corroborated 



291. 



Due to the importance of SFD FPTD in biological (and other) systems, there is a sur- 
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prising lack of research into the exphcit nature of the FPTD for a tracer particle of a SFD 
system, but a few exceptions exist, for example Li et al. used the Method of Images 
30| to help explain the FPT of proteins on "road-blocked" DNA strands in conjunction with 
3D diffusion. Taloni et. al. jsij, performed homogeneous SFD simulations, with particular 
emphasis on the long-time asymptotics. 



In t 
fusion 



le field of stochastic processes there exists a phenomenon known as anomalous dif- 
32M35I whose ensemble-averaged mean-squared displacement is represented by non- 
linear time dependence: (x^(t)) ~ t^^, where 2H 7^ 1. A prominent example of anomalous 



diffusion is the process of fractional Brownian motion (fBm) [36|, |37|. fBm is defined by a 
single parameter, the Hurst exponent H G (0, 1); and is a generalized, zero mean, stationary, 
Gaussian process with increments obtained from fractional Gaussian noise (fGn), produc- 



ing a position auto-correlation function |38'|: (i?//(t)-B//(s)) = C {t^^ + .s^^ — \ t — s P^), 
where C is a constant. In terms of displacement, a fBm time series is represented by 
the cumulative sum of the fGu, ,(«). whc^e^auto-eonelafou fuucfou is („(*),(0)) = 
H {2H - 1) + 2{l + 2H) 5{tf~'^", see [38]; where it is interesting to note that this 

or t > is equivalent to the memory kernel of the fractional Langevin Equation 



function 
(fLE) 
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40|. When H < 1/2, the pre-factor in the fGn auto-correlation function be 



comes negative for t > 0, illustrating an anti-correlated process 



in fields such as: Hydrology 
its name), diffusion o: 
genetic algorithms 



4l| . fBm is of relevance 



42| (the investigation from which the Hurst exponent derives 



biomolecules j43|, quantitative finance 4^ and even considered in 
45| . It should be noted that there are other processes which can give 
rise to anomalous-type diffusion, two of which are covered extensively within the literature: 
Continuous-time random walks 46|; and diffusion on fractal s | 47l| . For a comprehensive 

331, and later, Metzler and 



overview of anomalous diffusion see the reports by Metzler 
Klafter IsSj. 



The full functional form for the FPTD of a fBm process has remained elusive, but recent 
work into various FPT scenarios has become common of late. Such as, escape problems 
with fGn m 



and FPT of fBm in 2D wedge domains |49| . In ID, it was proposed by Ding 
aud Yaug M, shortly afte, a simple, sucemet. phys.eal argumeut was proposed by Kn,g 

and co-authors [511, and later rigorously proven by Molchan 
for the FPTD is 



fit) ^ t 



H-2 



52|, that the long-time form 

(1) 
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This result is a corner stone of our investigation. 

It is well-known in the mathematics literature that tracer particle dynamics in a hard-core 
atti ce g as (symmetric exclusion process) is equivalent to fBm, see for instance 22| (see Refs. 



16 



3l| for the case of continuum motion). Thus these processes/FPTDs are intimately 



linked via the theoretical long-time benchmark that is Molchan's result, Eq. ([T]). It was 
demonstrated in Ref. 16| that a homogeneous SFD system can be modeled with a Hurst 
parameter equaling i7 = 1/4, in the asymptotic time limit. We here go beyond previous 
studies, using Eq. ([T]) as guidance, by examining the full functional form of the FPTD, 
for SFD (homogeneous and heterogeneous) and fBm systems for all times, by extensive 
simulations and analytic approximations; so as to attempt to elucidate the entire FPTD 
relationship between these systems and the universality class to which they both belong. 

The layout of this article is as follows: In § [TTl the SFD tracer particle/fBm probability 
density function (PDF) is presented, and using this result, two Markovian techniques - 
Method of Images (MIA) and Willemski-Fixman (WFA) - are used to approximate the FPTD 
of the tracer particle within homogeneous and heterogeneous SFD systems, including also the 
FPTD of a fBm particle. The results are analyzed and compared to our simulations in the 
penultimate section. Whereas the MIA was used before, our explicit analytical inverse of the 
so called Willemski-Fixman relation is to our knowledge, new. Within this section, a simple 
functional form is conjectured to model the FPTD for all times and for all H. § 11111 describes 
how the SFD (homogeneous and heterogeneous) and fBm systems are simulated and what 
sets of parameters are used. Previous theoretical results for tracer particle dynamics are 
briefly reviewed. In § llVl the results concerning all aspects of the investigation are discussed 
in detail; this leads to § |V] that concludes the results herein, with a brief discussion on 
proposed further future investigations based upon this work. The more technical details 
regarding the derivations of the FPTD for both approximations are left for the Appendices 
lAl and [HI Data fitting procedures are situated in Appendix O 



II. APPROXIMATING THE FPTD 



To investigate the FPTD of tracer particle dynamics in SFD and fBm systems, which 
are inherently non-Markovian, we begin by analyzing the applicability of two well known 
Markovian methods; both of which require the following results: The tracer particle PDF 
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for SFD systems and the fBm particle PDF is Gaussian 12| (for infinite systems): 

1 / (x — Xo)^ 



P{x,t\xo) 



where Xq is the position of the tracer particle at the initial time, and x the position. The 
mean-squared displacement (MSD) is denoted S(t) and is, in the long-time for tracer particle 



(2) 



dynamics in SFD systems and of a fBm particle 



S{t) = {[x{t)-xof) = 2Ct 



2H 



(3) 



The angled brackets denote the ensemble average [53], and C is the effective diffusion con- 
stant. The explicit expression for C for SFD systems is discussed in § IIIII For fBm, C is 
the amplitude of the process (here as position) auto-correlation function, although often 
set to C = 1/2 in the literature 38|]. With these two expressions we can now develop two 
approximations for the FPTD. 



The Willemski-Fixman Approximation 

In a Markovian system, one is able to acquire the PDF and in turn, the FPTD through 



the well-known Willemski-Fixman (also known as Renewal Theory) method [5^. Simply 
put, a particle must have a first passage through point Xc at some time t' to reach point x 
at time t > t', i.e., functionally as 

P{x,t\xo)= [ f{x,,t'\xo)P{x,t-t'\xc)dt', (4) 

which constitutes a convenient integral equation for the FPTD, f{xc, t\xo), if one sets x = Xc 
in Eq. (jlj). However, due to the power-law argument with respect to time in the exponent 
in Eq. (jlj) for the PDF^Eq. ([2]), we cannot use the usual Laplace-transform ;echnique 
for inverting Eq. [54]. However, through the use of Mellin transforms 55|, we are 
able to follow the usual line of attack in this type of problem: transform the convolution 
into a product in "frequency" (Mellin) space p, rearrange to find the Mellin-transformed 
FPTD f{xc,p\xQ), then transform back to t— space and re-acquire f{xc,t\xo). The explicit 
derivations are left to Appendix |3 and the final results are given here: 

Using a series expansion approach [56] of our exact expression for f{xc,p\xo), we find 
an exact expression for the FPTD within the WFA for the general case (0 < if < 1; with 



H 7^ H^^, see below) through the inversion of the Melhn transform. Exphcitly, 



a 



-1/{2H) 



T{l-H) 



n=0 



n (-ir r(^) 

n\ r (i7 - 1 - n) 



+ 



m J2 mt)] 



2Hm {-irT{l-H-2Hm) 



m=l 



ml 



r {-2Hm) 



(5) 



where a = (Ax)^/4C, with Ax = Xc — Xq, T{z) is the Gamma function 57|, and K{t) is 
defined as 



a 



1/2H 



t 



(6) 



Eq. ([5]) is very computationally efficient for calculating FPTDs, but is only valid when H 
is not a rational number of the form H*^^ = (n + l)/(2m + 1) where n and m are positive 
integers (see Appendix |X] for further details). Eq. generalizes the long-time expression 



m 



58] to include also shorter time-scales. Eq. is equivalent to the usual Brownian 



motion result when H = 1/2; again see Appendix [Aj Alternatively, the inversion back from 



Mellin-space can also be completed using Weyl fractional derivatives 
an equivalent expression 

2sin(7riJ) 



32l . l55| . which yields 



fwFiXc,t\xo) 



n 



(7) 



POO 

I{z)= / y''[2H{z + yy'' + l-2H]iz + yf''-'exp[-iz + y)''']dy, (8) 
Jo 

vahd except for K{t) -> for < 1/3. 

To check the validity of the WFA, we look toward the long-time limit, t ^ a^/^'^^\ for 
comparison to the theoretical result given by Eq. ([T]). In this limit where K{t) — )■ 0, the 
series expansion result gives: For i7 > 1/3 



fwF{Xc,t\xo) 



{2-H) ^ 



(9) 



For H < 1/3 



fv^F{Xc,t\xo) 



2HT (1 - m) 



V{l-H)V{l- 2H) 
When H = 1/3, the long-time limit becomes 



a 



-1/{2H) 



[Kit)] 



(10) 



f^NF{Xc,t\xQ) 



(11) 
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as shown in Appendix[Xl Euler's constant 57| is given as 7^;. Note that our long-time results 



coincide with Bologna et al. [58]; but therein, the authors solve only up to a pre-factor, 
whereas our result provides a general explicit solution including the pre-factor. Combining 
Eqs. (Q, (do]), and gives the FPTD within the WFA for all H in the asymptotic limit. 
Interestingly, when H > 1/3, as in Eqs. Q and f lTT]) the WFA approach is in agreement 
with Molchan's result. However, for the case that H < 1/3, Eq. (ITUI) shows a different 
exponent to Molchan's {H — 2). Thus our work shows that i7 < 1/3 is, in a sense, "too far" 
from the Brownian case {H =1/2) for renewal theory to apply. Further analysis of these 
results are left to § [TVi 

The Method of Images Approximation 

The Method of Images is a technique used in calculating properties of electrostatics jsO], 
but has also been applied to Markovian stochastic systems [1]. This technique calculates 
the resultant PDF from having a PDF in the absence of a boundary, centered at xq, and an 
identical PDF centered at 2xc — Xq, which is inverted. The result "creates" the existence of 
a boundary at x^, from this, the FPTD can be calculated. For the case of non-Markovian 
systems, such as fBm or the SFD system, this method is only an approximation, due to the 
fact that the state of the system is dependent upon all previous trajectories; namely the 
system is not memory-less. Never-the-less, this technique has been used in the literature 
to describe the FPTD in similar situations 29||. The resultant FPTD (see Appendix [B] for 



derivation) is 

/Mi(a:c, tW) = ^7^^^ exp [-^^ j ■ (12) 

The FPTD for a Brownian particle is produced when H = ^ and C = D, where D is the 
diffusion constant, as it should. However, for any system such that H ^ 1/2, the long-time 
behavior predicted by the MIA, /mi(^) ~ t~'^^~^^\ where this result and Eq. ([T]) are in clear 
conflict. A more in-depth analysis of the MIA with respect to simulations is left to § IIVI 

Conjecture for the FPTD 

As we will see in § IIVI the Markovian approximations fail to produce the correct full 
functional form FPTD for the two physical SFD systems present here, and more explicitly, 
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analytically, only the WFA appears to capture the theoretical long-time result for H > 1/3. 
Therefore we here provide a simple conjecture made in light of the results in previous work 
and more simplistic systems to account for the full form. 

Our proposed function for the FPTD has two main features: (i) a function which becomes 
exact for the Brownian case and (ii) approaches the correct long-time behavior for both SFD 
and fBm systems, Eq. ([T]), for arbitrary H (which, obviously, includes those values which 
represent SFD homogeneous and heterogeneous systems). Namely 



/c(xc,t|a;o) = Vtt" ^exp I -7 



with dimensionless fitting parameters 7 and (5. The normalization constant, is 



(13) 



l-H 



2C 



7 



l-H 

2HI3 



2HI3 

Consideration of this conjecture in conjunction with simulations is left to § IIVI 

III. SIMULATIONS 

Simulating SFD: To computationally simulate the FPTD in a single-file system, the 



Gillespie-type algorithm for hard-core lattice dynamics presented in Ref. [2l| is implemented. 
Stochastic time series are generated for both homogeneous and heterogeneous population 
types through the following steps: (1) Place the particles in their initial, thermally equili- 
brated, positions (the tracer particle is positioned in the middle of the system with equally 
many particles randomly distributed to the left and the right). (2) Move a random particle 
according to the algorithm in (2l| and update the time t. (3) repeat step (2) until either 
t > 4top (some designated stop time) or the tracer particle position becomes > Xc, and 
record the corresponding first-passage time; see Fig. [T] for a diagrammatic explanation. 

Both homogeneous and heterogeneous particle populations were simulated using the pa- 
rameters provided in Table [H The simulation box size, L, was chosen such that the results 
accurately displa y th e dynamics for a semi-infinite system given the finite nature of a com- 
puter simulation |59|. The lattice simulations conducted have focused on one particle density 
(p = 0.25), one particle size, a (which is equivalent to the lattice site length), and four dif- 
ferent absorption point differences. Ax (shown in Table [IT]) . 
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FIG. 1: Schematic representation of the lattice simulations of a SFD system. All particles (including 
the tracer - here in black) move under Brownian motion and are hard-core (the particles cannot 
occupy the same lattice site), meaning that they cannot pass each other, keeping their order for 
all time. Hence the tracer is in the center of all other particles for all times. The top panel shows 
the start of the simulation in thermal equilibrium, whereas the bottom panel shows that after some 
time, the tracer has achieved a first passage event. 



TABLE I: SFD Simulation Parameters. Not applicable is abbreviated to N.A. 



Parameter 


Homogeneous 


Heterogeneous 




systems 


systems 


System Size, L (lattice site width a) 


10^ a 


2 X 10^ a 


Particle Number Density, p 


0.25 


0.25 


Diffusion Constant, D 


1 a^s 


N.A. 


Average Diffusion Constant, D 


N.A. 


1 a^s 


Stop time, t > tgtop 


2.5 X 10'^ aV-D 


10^ a^/D 


Friction constant PDF exponent, a 


N.A. 


0.5 


Initial tracer particle position, xq 








Simulation Ensemble size, N 


2.5 X 10^ 


2.5 X 10^ 



For heterogeneous SFD systems the particles' friction constants (^ oc D ^ in units of 



thermal energy; D is the diffusion constant) are drawn 



rom a so-called heavy-tail (HT) 



power-law distribution (the first moment is not finite) jl8l. |60|. The HT frictional constant 
distribution is chosen as g{C,) = S,^^^^", for large C,, where < a < 1 and H is a normalization 



constant. For such heterogeneous SFD systems the dynamics are found to be inherently slow, 
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and therefore large systems and simulation stop times are required (compared to the same 
system with a homogeneous population), this is accounted for in our study, as shown in 
Table E 

Simulating fBm: We here use the so-called Spectral simulation method; as it uses 
the spectral density of the fGn increments to calculate the time series. This method was 
chosen for its speed (use of fast Fourier transforms) and simplicity of implementation, which 
approximates well fBm, for long series |6l|. For sub-diffusive systems with slow dynamics, a 
relatively large (with respect to SFD systems herein) effective diffusion constant, C = 5, is 
used to efficiently calculate and display all relevant time frames in a reasonable computation 
time. See Table UTTl for all fBm simulation parameters j62|, and results. 

Collapsing Data: When comparing sets of data from the same system (characterized by 
H) with different absorption points and/or generalized diffusion constants, it is convenient 
to collapse the data to remove these dependencies on Xc and C and thereby be able to 
extract a universal functional form for the FPTD. This is done through knowledge of length 
versus time scaling in anomalous diffusion, provided by Eq. ([3]). Through observation of this 
equation it is convenient to introduce a scaling factor w = (Ax)^^^^ C^^'^^ with dimension 
inverse time, and to scale FPTDs with this factor, as it is done in Figs. [2]to|5l 

For homogeneous SFD systems it has been shown (l?! that the MSD for a tracer particle 
is described by, see Eq. Q, 

H.\. (14) 

p V TT 4 

where D is the single particle diffusion constant, and a is the size of the particle. For 



heterogeneous SFD systems Lomholt et al. 
factor and Hurst exponent to be 



18] find the approximate long-time MSD-pre- 



C-^^^, H=^^,, (15) 
2^T{1 + Sy 2(l + «)' ^ ^ 



where k = p'^ksTil - pa)-'^, x = {4:k)^'^^{E7i/ sm[2HTr/{l - 2H)]Y'^^-^^\ and a is the 
exponent in the HT friction PDF. Note that the Hurst exponent satisfies H < 1/4 in 
these heterogeneous systems, as < a < 1; we use a = 1/2 during this investigation (see 
Table [T]). From Eq. (fT5|) . to may be calculated for heterogeneous SFD systems, likewise for 
homogeneous systems using Eq. ([1] 
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IV. RESULTS 



SFD simulations: Each different data set (defined by its unique absorption point, Xc) 
is sliown in a collapsed plot in Fig. |2]for homogeneous particles and Fig. [3] for heterogeneous 
systems. There are concurrent motifs in both FPTD plots; namely the short-time concave 
characteristic maximum moving into a long-term power-law type structure. For the hetero- 
geneous case, Fig. [21 there is more variance in the short-term time regime compared to the 
same time regime for the homogeneous. Fig. [21 This is because the system is given a new 
.andom set ot diffusion constants fot every singulation (the so-called heterogene.ty-avetaged 
case [18]), therefore giving a super-position of non-averaged FPTDs (see discussion at the 
end of this section) as demonstrated in Fig. [3l inset. 

Long-time SFD asymptotics: An important check of simulating homogeneous SFD 
data is that it should agree with Eq. ([11), for long times. Upon evaluation of Fig. [21 where 
Molchan's asymptotic result is placed to guide the eye, one sees that the data agrees well 
with theory, in accordance with previous literature findings [31] for such a system, further 
confirming that SFD systems and fBm belong both to the same universality class in the 
asymptotic limit. 

A novel result then was to verify if Molchan's work agreed with the heterogeneous SFD 
system proposed and investigated by [18], assuming H is given by Eq. f[T51) . In Fig. [21 in 
both the main and inset plots, the theoretical asymptotic result (which is again placed to 
guide the eye) agrees well with the long-time FPTD for heterogeneity-averaged simulations. 



SFD versus fBm: In both collapsed plots (Figs. [II'G^j and [31 [65]) the fBm simulations 
agree well with the SFD data, including the scaled power-law pre-factor. 

Comparison of the short-time regimes in Fig. [21 and in particular the inset, illustrates 
that the fBm and homogeneous system have the same FPTD on all temporal scales. 

Within the short-time heavy-tailed heterogeneous SFD regime, although possessing the 
same structure, the fBm simulations do not coincide with the heterogeneous data (see dis- 
cussion at the end of the section). 

SFD data versus approximations: In the homogeneous system, as demonstrated in 
§ ini the long-time MIA FPTD, Eq. ([T2D, does not agree with Eq. ([11). With inspection of 
Fig. [21 we note that although the concave structure of /mi(^) is present and similar to the 
data, these two FPTDs do not coincide, making for an obvious case that fui{t) does not 
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FIG. 2: Collapsed (log-log) plot of FPTD for a homogeneous SFD system with different absorption 
points and fBm FPTD (for simulation parameters see Tables HI and IIIH respectively). Immediately, 
the collapsed plot shows that homogeneous SFD and fBm have the same FPTD dynamics over all 
time frames in the correct scaling, ^ IIIII The MIA, Eq. (|12|) . and the WFA, Eq. (0), are shown for 
comparison, both of which show poor agreement. The averaged proposed functional form, /c(t), 
Eq. ()13p , is constructed by collapsing all simulated data and fitting (see Appendix [C]) , keeping 
the power-law exponent in the pre-factor fixed to H — 2; with H = 1/4, and setting C according 
to Eq ()14p . The parameters (Table |TTj) were then averaged and a single mean curve plotted. This 
conjecture shows excellent agreement with both anomalous diffusive systems on all time scales 



n 

63|. The Molchan long-time prediction [52] is given to guide the eye. The fBm FPTD consists 



6j| and displays 



of two different absorption points {Ax = 50, Ax = 100, 6 x 10^ simulations) 
good agreement with SFD results. INSET: The short-time regime (linear axes) agreement between 
homogeneous SFD, fBm, and fdt). Remaining simulation details are presented in Tables HI IIIH 
andini The subscript s on each axis variable denotes "scaled", namely tg = wt and fs = w^^f{t), 
where ro = Ci/2^ (Ax)-^^" . 
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FIG. 3: Collapsed plot of FPTD for a heterogeneous SFD system (Table[Ij) with different absorption 
points. Also within, fBm FPTD is plotted for H = 1/6, see Table UTll The long-time dynamics for 
both systems agree very well with each other and also with Molchan's equation, Eq. ([1]). In the very 
short-time the systems part, most likely because of the complex nature of heterogeneity- averaged 
SFD systems (see inset). The MIA, Eq. (fT2|l . and the WFA, Eq. ([5]), are shown for comparison. 
Both approximations show ill agreement with the simulated data for both systems and the theory. 
/c(t), Eq. p3|) . was fitted to the fBm data (displaying excellent agreement). Tabl e IIIH as opposed 
to the SFD data, due to its poor fit (discussion in inset caption here, and § HYj) j65||. The data is 
scaled (signified by subscript s) using Eq. (fT5|) . with H = 1/6 (since a = 1/2) - see Fig. [2] caption 
for further explanation of the scaling. All remaining simulation details are presented in Tables [I] and 
mil INSET: Collapsed plot of the FPTD for 3 different sets of heterogeneous friction constants, kept 
constant for each simulation. For this non-averaged case we use Xc = 50 for all simulations, with all 
other system parameters displayed in Table |TJ The inset illustrates the fact that no self-averaging 
takes place in this heterogeneity-averaged system (see § IIVI for further discussion) . 

agree with the SFD homogeneous system on any temporal scale. Similarly for the MIA in 
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Fig. El 

Concerning both Figs. [2] and El the WFA FPTD also produces the required short-time 
concave structure and hes closer to the simulated data than the MIA (remembering the log- 
log scale), but still does not agree with it. In the long-time, as visible through Eq. (ITOl) and 
Fig. [21 we see that although the power-law structure of /wf(^) does have a closer gradient 
than the /mi(^), we can still conclude that clearly the fwF{t) does not, again, agree with the 
simulated data on this scale. 

fBm and WFA: When one considers sub-diffusive fBm, under the criterion: 1/3 < H < 
1/2, we see, via Fig. H] (where H = 0.35, top, H = 0.45, bottom), that the WFA agrees well 
with the fBm data: in the early time regime, and in the long time, has the correct gradient 
and a power-law pre-factor becomes exact as — )■ 1/2 (Brownian motion), as it should (and 
so to does the MIA in this limit). 

For super- diffusive fBm, H > 1/2, our simulated results, with H = 0.75, except pre- 
diction of the correct long-time power-law exponent, show poor agreement with the WFA 
approximation, on every scale - see Fig. |5l For simulation and fitting data see Table IIIII 

Conjecture results: The viability of the trial function, Eq. (fT3l) . was measured with 
comparison to simulated data (both SFD and fBm; Tables HTl and [XTTt respectively). 

In the homogeneous system. Fig. [21 the power-law exponent, is fixed to the theoretical 
result: {H — 2) ; where H is taken from the literature in the form of Eq. (1141) . leaving only 
two parameters, 7 and /3; to be fitted (see Appendix [C|), resulting in 

Table [III H- 

The low 

TABLE II: Homogeneous fc{i) fit parameters, is the normalized Chi-squared parameter. Raw 
data placed into 30 natural log-bins (see Appendix [C]) . 





7 


/3 




25 


0.72 ±0.03 


0.96 ±0.02 


1.49 


40 


0.67 ± 0.03 


0.97 ±0.02 


1.64 


50 


0.64 ± 0.03 


0.98 ± 0.02 


1.61 


75 


0.60 ± 0.02 


1.00 ±0.02 


1.33 



value, shown in Table [TTl indicates that our proposed form for the FPTD is indeed an 
accurate one. Further, the results in the /3 column of Table [Til leads us to propose that /3 = 1 
(as is for the Brownian case), i.e. the FPTD is well approximated by the propagator, Eq. 
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10"2 lO^'' 10° 10^ 10^ 10^ 

ts 



FIG. 4: Collapsed plot of FPTD of sub-diffusive fBm with H = 0.35 (top), H = 0.45 (bottom) 
(simulation parameters: Table IIIip . Both panels show that the WFA has the correct heavy-tail 
gradient (compare to Molchan's prediction). The fBm also agrees with Eq. ([T|), as the fdt) tail is 
fixed (H — 2); and the data is modeled well by our conjecture, Eq. (|13p . Our conjecture models all 
time scales well, see Table UTTl for quantitative details. 

([2]), multiplied by a power-law pre-factor chosen to reproduce Molchan's long-time relation. 

For the heterogeneous case, the simulation data have been subjected to the same treat- 
ment as applied to the homogeneous system data; yielding poor results. The full func- 
tional form of the heterogeneous system is more complicated and not amenable to expla- 
nation by our simple conjecture. This complexity, as stated previously, arises predom- 
inantly through creating an entirely new heterogeneous population for each simulation 
(heterogeneity-averaged). Each different population, three of which are illustrated in the 
inset of Fig. |3l has a different short-time FPTD structure relative to other populations as 
well as a different pre-factor in the power-law tail. Thus, in short, the FPTD for heteroge- 
neous SFD shows the same kind of lack of self-averaging as displayed in the MSD for this 
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FIG. 5: Collapsed plot of FPTD for super-diffusive fBm, H = 3/4, see Table |TTI1 Both the WFA 
and the fBm data only agree with Eq. ([TJ in the long time, as expected. Using the conjecture, Eq. 
(jl3p . keeping H fixed, the 2 remaining degrees of freedom (7, /3, see Table UTT)) cannot account for 
the discrepancy seen between the our conjectured FPTD and the fBm data. It is apparent that the 
MIA fails on all time scales. INSET: Collapsed plot on linear axes illustrate short-time dynamics 
(fBm data as crosses). 



type of system [18|. In particular, the pre-factor in the long-time asymptotics is a random 
quantity dependent on the particular realization of the friction constant population. 

Upon inspection of Figs. |2] to SI it is easily seen that the conjecture, Eq. f|T3|) . works 
well at modeling the FPTD dynamics of fBm for sub-diffusive and diffusive systems. When 
considering super-diffusive motion. Fig. [5l it is clear that the conjecture is not adequate to 
deal with these dynamics. Quantitatively, these statements are displayed in Table IIIII 
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TABLE III: fBm Simulation and Fit Parameters. Each simulation has: C = 5; tgtop = 10''; 

= 6 X 10^. Raw data placed in 50 natural log-bins, fitted to parameters (7,/3) (see Appendix 
IC]) to test the applicability of our simple conjecture, Eq. ()13p . is the normalized Chi-squared 
parameter. 



H 


Ax 




/3 


y2 

A 


1/6 


25 


0.726 ± 0.004 


1.00 ±0.02 


1.12 


1/6 


75 


0.556 ± 0.004 


1.004 ±0.002 


1.21 


1/4 


50 


0.639 ± 0.005 


0.987 ±0.003 


1.70 


1/4 


100 


0.595 ± 0.005 


0.984 ± 0.003 


1.65 


7/20 


50 


0.667 ±0.006 


0.958 ± 0.005 


3.16 


9/20 


50 


0.595 ± 0.007 


0.973 ± 0.006 


0.985 


3/4 


200 


0.0862 ± 0.0004 


1.70 ±0.03 


55.7 



V. SUMMARY AND DISCUSSION 

We investigated the first passage time densities (FPTD) of a tracer particle in a single- 
file system with two different population types; homogeneous (all having the same diffusion 
constant), and heterogeneous (friction constants drawn from a heavy-tail power-law distri- 
bution). Along side, the fractional Brownian motion (fBm) FPTD was investigated. 

Theoretically, two methods were used to approximate the full functional form for FPTD 
analytically: the Method of Images approximation, Eq. (I12p . and the Willemski-Fixman 
approximation, Eq. Moreover, a conjectured form for the FPTD, Eq. (1131) . was in- 

troduced. Numerically, the SFD and fBm systems were simulated stochastically with the 
Gillespie-type algorithm presented in Ref. j2r|, and spectral approach joil, respectively. 

Our main conclusions are: 

• The MIA derived here, and previously used in the literature, does not approximate 
the FPTD at any temporal scale. 

• With the use of Mellin transforms we have found an exact result for the full FPTD 
within the Willemski-Fixman approximation. To convert the inverse FPTD from 
Mellin frequency space two methods were used: the series expansion approach and 
the Weyl fractional derivative approach. The WFA does not agree with simulated 
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the SFD FPTD (homogeneous and heterogeneous) nor the fBm FPTD for H < 1/3. 
For sub-diffusive and Brownian motion, i.e. 1/3 < H < 1/2, the WFA approximates 
the FPTD well at all times, including the theoretical long-time slope, Eq. ([1]), (but 
does not predict the correct pre-factor) becoming exact when H = 1/2. In the super- 
diffusive regime, H G (1/2, 1), the WFA manages to capture only the long-time power- 
law exponent correctly. 

• We show through simulations that the FPTD for the homogeneous SFD system is 
equivalent to the fBm FPTD with H = 1/4, for all times, in the correct scaling. 

• We find that the FPTD for a heavy-tailed heterogeneous SFD system and the corre- 
sponding fBm share the same power-law exponent for long times (given the correct scal- 
ing, Eq. (fT5|) ). The general heterogeneous SFD FPTDs show a lack of self-averaging. 

• A simple conjecture, Eq. (IT^ . is proposed and adequately show to model the full 
functional form of the fBm FPTD for H < 1/2; and also that of the homogeneous SFD 
FPTD. Eq. (1131) only captures the asymptotic power-law exponent for heterogeneous 
SFD systems and super- diffusive systems where H > 1/2. 

Our work presented herein also pertains to research on FPTs in many-body "swarm" sys- 



They have investigated the 



tems; such as the recent work by Mejia-Monasterio et al. 
FPT for a search by a swarm of independent searchers. Our results initialize the general- 
ization of this investigation, and further paves the way for the understanding of FPTs in 
general interaction, non-independent, many-body swarm systems. 

Whereas there is a wealth of knowledge about first passage problems for Markovian 
systems |l|, our understanding of the corresponding problem for non-Markovian dynamics 
is far from complete. This study shows quantitatively the limitations of two commonly used 
approximations, provides extensive simulation results for homogeneous and heterogeneous 
system, and a conjectured form for the FPTD for the fBm, thereby providing headway into 
future studies. 
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Appendix A: WFA Derivation 

The canonical method for solving the Willemski-Fixman problem, Eq. (jl]), is to use 
the Laplace transform (operator L) to take advantage of the convolution property of this 
transform. However, this transform for the function: L [exp {—Vt~^)] with respect to t; 
does not exist in closed form for arbitrary r. In here, we instead therefore use the Mellin 
transform technique. 

Beginning from the integral Eq. (j3]), plugging in the propagator, Eq. ([2]); multiplying 
both sides by , we write 

exp = / Fix,,t'\xo)Git/t') — , (Al) 

Jo t 

with F{xc,t'\xQ) = t' f{xc,t'\xo), where a = (Ax^)/(4C) as before, and 

where B(t) is the Heaviside function. We now can solve Eq. (lAip . using the Mellin Trans- 
formation (see 55| for more details). This transform is defined as 

POO 

gip) = M[git)]^ / t^-'git)dt. (A2) 

^0 



Using the following identities: 



M 



M[9it-'")L>o = b(^]^ (A3) 



1 - / ~p ' 



w \ w 



the formal solution to Eq. flAip becomes 

^^a^^'^'T (^^^ = F{x^,p\xo)G{p), (A4) 



2H 

where 



G{p) 



/OO POO 
tP~h^{t - l)-^Q{t - l)dt = j tP-^"-\t - l)-"dt = B{-p, 1 - H). 
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B(z,w) is the Beta function 



57| . where B{z,w) = [T{w)T{z)] /T{z + w). Following this 



route, one finds that Eq. (lA4p becomes 



F{xc,p\xo) 



2HT{1 - H) T{-p) 



(A5) 



for Re(p) < 68|. Using the fact that M [F{t)t = F{p — 1), (remembering F{xc,t\xo 
t/(xc,t|xo)) one finally gets 

,{p-i)/2H Y{2-H-p)T{^) 



f{Xc,p\xo) 



(A6) 



2HT{1-H) T{l-p) 

for Ke{p) < 1 , which gives the closed form expression for the Mellin transform of the first 
passage time density f{xc,t\xQ), within the WFA. 



1. Series Expansion Approach 

From Eq. (1A6|) . and the formal Mellin- inversion formula, we have: 

f{xc,p\xo) = -^^——^— / rVwr(2-/J-p)— ^t/p. (A7) 
2HT (1 - H) 2m J^_i^ T (1 - 

We can deform the original contour where c is chosen such that the integration is over the 
fundamental strip (here Re(p) < 1), so that it becomes a square, side length R, which 
encompasses the fundamental strip, encloses all poles in the domain Re(p) > 1 and has a 
negative orientation (anti-clockwise) - see [56], Appendix A therein. We can then use the 
residue theorem to compute Eq. (1A7|) : 



2HT (1 - H) 

Residues 



t-v^r(2-i/-p)^A^ 



(A8) 



where we used the fact that the integral along the upper, lower, and right parts of the square 
vanish as R ^ oo. Note the minus sign in the pre-factor due to anti-clockwise contour. T{p) 
has simple poles at p = 0,-1,-2,...; and the function inside the square brackets on the 
RHS of Eq. flA8|) therefore has (potential) poles at 

Pn = 2-H + n- n = 0,l,2,... (A9) 

Pm = l + 2Hm; m = 0,l,2,... (AlO) 
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To proceed we need to know the behavior of F 
for a simple pole the residue is calculated as: 



close to Zm = — m/3. We have that 



Res 



a — z 



(3 



lim (Z — Zrr,)T 



a — z 



(3 



after some manipulation, by inserting Zm = (x + and using T{z) = z ^r{z + 1), we find 
that 

-P -P -13 , 



Res 



r 



a — z 



/3 



m/3 (m-l)/3(m-2)/3" 



-/3)r(i) 



ml 



-a. 



(All) 



Noting that \jY{z) has no poles, and combining Eqs. flASP and (lAlip we get 



-1 



2EV (l-H) _ 



n=0 



+ 2_^t ^"cr 2H A L — 2H- ^ 



m=l 



r(i-p. 



substituting in K{t), see Eq. ([6j), and the values for p„ and p^, one arrives at the exact 
expression for the first passage time within the WFA, which can be used to numerically 
evaluate the FPTD: 



a 



-1/2H 



V{l-H) 



2H 



:-!)" r(^) 

n\ V{H -l-n) 



^K(f) V \K(f)]'^- {-ir T{l-H-2Hm) 
+m ^ [Kit]] --^ r (-2Hm) 



n=0 



(A12) 



T(-2Hm) 

We point out that we above assumed simple poles, i.e. that p„ ^ pm- Using Eqs. flAOp and 
(1A10|) . we see that we therefore require H ^ {n + 1)/ (2m + 1) for Eq. (1A12() to be valid. 

Now, an immediate check for the calculation of Eq. (IA12P is performed when H = 1/2, 
i.e. the system's dynamics is Brownian motion. Then 



f{Xc,t\Xo) 



a 



-11 °o (—AY 

Try [mf E [^(^)]" + E [^w] 



-i)™r (i -m) 

m! r i—m) 



r(i) , 

V2/ \ n=0 m=l 

where we note the sum of n term is the power series of an exponential function, and that, in 
the second sum, over m, l/r(— m) = for an integer m > 0. Substituting in the functional 
value for K{t), we get 

rt '^/^exp 

as required [l|. 



TT 



-a 
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2. Asymptotic Time Limit: K{t) — ;> 



In the long-time limit n = term and the m = 1 term in Eq. (1A12P give the dominant 
contribution to the FPTD, namely 

We see that the exponents in Eq. (1A13P are equal if if = 1/3. Thus this obliges us to 
make three distinct cases for the FPTD: when H > 1/3 we use the reflection formula for 
F— functions and arrive at Eq. ([H]); in the case that H < 1/3, we flnd Eq. (ITU]) , using 
F(l + z) = zr{z). For the case that H = 1/3 (or more generally of the form if*^ = 
{n + l)/(2m + 1)), a double pole exists and certain care is needed. 

Double Pole: When H*^ = (n + l)/(2m + l), e.g. H = 1/3, one or several double poles 
exist and should be dealt with accordingly. To see so, if we substitute if = 1/3 into both 
Eq. and Eq. (llOp we flnd a diverging pre-factor which is proportional to 1/ {H — 1/3), 
which is not correct as f{xc,t\xQ) is bounded, namely: f{xc,t\xQ)dt = 1. In order to 
resolve the double pole issue we need to revisit Eq. (lASp . 

The residue for a pole of order 2 is computed as 

Res [F{z)] = hm [{z - z^YF{z)] . (AM) 

z^zm az 

Considering Eq. (lASp . we note that we need to evaluate - following the same steps that led 
to Eq. (lASp - (remembering: Zn = f] + n, = 0, 1, 2, ... and = a + m/3). 



Res 



T{ri-z)T 



a 



(3 



a—z-\- 



X 



{z — a){z — a — j3) ■ ■ ■ {z — a — {m — l)/3) 
T{ri- z + {n + l)) 



(A15) 



{z — ri){z — 7] — 1) ■ ■ ■ {z — 7] — {n — 1)) 

Eqs. (1A8I) and (lAlSp allow us to compute f{xc,t\xo) (in principle), also for n and m values 
satisfying H*^ = {n + l)/(2m + 1). For a given (n,m) from Eq. (lAlSp . we see that the 
residue for the double pole involves 2 + n + m terms. Let us now limit our derivation to the 
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analysis of if = 1/3 and long times (where n = and m = 1), then Eq. (lAlSP becomes 



Res 



T{r]-z)T 



a — z 



/3 



X 



OL — Z 

/3 



+ 2 r(r7-2 + l) 



z — a. /3 



+ 



z — a 

V'{ti-z^\) 
Y{j]-z^\) 



For the investigation in this report, we have /3 = 2if; 5^ = Pm = 1 + IHm — )■ 1 + 2if; and 
a = 1; therefore z„ = p„ = 2 — if, and finally t] = 2 — H . Substituting these values in above, 
and setting if = 1/3 we find 



Res 



r(2-ii"-z)r 



1 - z 
2H 



, 5 
1 + 37. 



where (r'(x)) / (r(x)) 



n. 



x=l 



jE (Euler's constant) [57]. Eq. flA8l) now becomes (with 



P = = 2 — H) Eq. (ITT]) , as given in the main text. 



3. Weyl Fractional Derivative Approach 



Let us now derive Eq. ([7]) in the main text. To obtain the FPTD, one has to invert the 
transform in Eq. (1A6|) . The inverse Mellin transform is defined as |55| 



m 



2iii 



f{p)t-Pdp, 



where a < c < b and (a, b) is the the fundamental strip; in this case, (—00, 1). Transforming 
the inverse integral of Eq. (1A6P to the variable p = 2 — H — pwe obtain: 



f{Xc,t\xo) 



T{l-H) 



[K{t)t''<j>[K{t)] 



2H 



r[p-{i-H) 



t 



(AI6) 

(A17) 
(A18) 



The quantity (j){z) can be written in terms of H- functions or Fox- functions, but from 55] , 
Eq. (8.5.23) one has 



T{p) 



r(p-/9) 



g{p-/3),z 



[g{z)] 



(-1)" 
r(n-/3)d^" 



{z' - zY~^-^g{z')dz\ (A19) 
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where is the Weyl fractional derivative and n is the smallest integer such that n — f3 > 



55|. In this case - Eq. (1A17P - one has /3 = 1 — H which implies 0</3<l, son = l and 



g{p) = r {p/2H). Therefore g{z) = M'^ [T {p/2H) , z] = 2H exp (--2^^), and 

(f){z) = W^-"" [exp (-2'^)] (A20) 
-1 d 



T{H) dz 
2H 



poo 

/ iz'-z)-^'-''^exp[-iz'f'']dz' 

J z 

POO 

/ y-('-^\z + yr''exp[-iz + yr]dy 
Jo 



T{H) 

where, in the last step, we changed variables y = z' — z. The diverging factor y~^^~^^ in 
the integrand is inconvenient for numerical evaluation of the integral, so to remedy this, we 
integrate by parts to get 

"^^'^ = fm 

Using the reflection formula for Gamma functions 57|, r(z)r(l — z) = it/ sin(7r2;), (note: 



POO 

/ y"[2Hiz + yy"-i2H-l)]iz + yy^-'exp[-iz + yy^]dy. (A21) 
Jo , — , 



< z < 1) we finally get - by inserting Eq. (lA2ip into Eq. (lA16p . Eq. ([7]) given in the main 
text. Careful considerations show that the integration by parts is only valid for K{t) 7^ 
when H < 1/3 providing limits for the validity of Eq. (jTj). 



Appendix B: MIA Derivation 



Given the result of the PDF from § |TT]- Eq. ([2]), we have the PDF in the presence of an 
absorption point located at Xc, 

Pmi{x, t\xo) = P{x, t\xo) - P{x, t\2Xc - Xo). 

This combination of PDFs enforce the absorbing boundary condition Pui{x = Xc,t|xo) = 0. 
Using this result, the Survival probability is Suiit) = S^it) — S^it), where 

(Bl) 



(B2) 



SAit) = / P{x,t\xo)dx; SB{t) = / P{x,t\2xc — Xo)dx. 
J —00 I — I J —00 

Making use of the Error function, Erf[z] |57|, one finds that 



Erf 



( 



Xc - Xo 



5f 



1 + 



Erf ( 



Xo - Xc 



The relationship between the FPTD and the Survival probability is /mi(^) = — ^5'Mi(i) jH, 
this leads to the final result predicted by the MIA for a non-Markovian system - Eq. fll2p - 
in the main text. 
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Appendix C: Data Reduction and Fitting 



The simulated data presented here was subjected to the following data reduction pipeline: 
first passage times are culminated from simulations then the data are log-binned and plotted 
as histograms. When fitting of histograms is required, the boot-strap method joQl is used to 
estimate the error of each bin. In this method one creates a large number (in our case 10^) 
"synthetic" raw data sets from the original raw data (the vector/set containing the FPTs). 
Each synthetic set is accomplished by selecting, at random (from a uniform distribution), 
an entry from the original set, with replacement, until the synthetic set being created has 
the same number of entries as the original. This is repeated until the desired amount of 
synthetic data sets is achieved. Each synthetic set is then binned into bins of the same 
size as the original data set, and the variance for each bin (acquired from binning a large 
number of synthetic data sets) is then used as the error for the original bin. This error is 
subsequently used to weight the least-squares fitting. 

After boot-strapping, the Levenberg-Marquardt regr ession (LMA) is employed to find the 



best fit 



70 



66(1, the LMA was employed to minimize 



■ l7l|. When fitting the FPTD conjecture 
the Chi-squared value, x^. To be more explicit, the normalized Chi-squared value, was 
minimized as it is more easily interpreted (namely without the use of tables). This is just: 
X^/dof, where dof is the degrees of freedom of a fit. If this value lies between 0.5 < < 1.5 
the model is assumed to be an accurate explanation for the functional form of the data. 
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